Data description

Data was simulated using VirtualCommunity code.
Simulated data contains 20 data sets.

NB - if chain/sample parameters is modified - should be changed in the whole section(dataset)

Environment filtering 5 species

JSDM

data<-sim_data$EnvEvenSp5
data <- list(
  Y = subset(data, select = -env),
  X = cbind(1, scale(poly(data$env, 2))),
  covx = cov(cbind(1, scale(poly(data$env, 2)))),
  K = 3,
  J = ncol(data) - 1,
  n = nrow(data),
  I = diag(ncol(data) - 1),
  df = ncol(data)
)

Y_cor<-cor(data$Y)

to_prec<-function(m){
  n<-dim(m)[1]
  Tau_n<-matrix(nrow=n, ncol=n)
  for (j in 1:n) {
   for (k in 1:n){
    Tau_n[j, k] <-  -m[j, k]/sqrt((m[j,j]*m[k,k]))
   }
  }
  return(Tau_n)
}

#Tau_n<-matrix(nrow=dim(model$mean$Tau)[1], ncol=dim(model$mean$Tau)[1])
Tau_n<-to_prec(me5$mean$Tau)
#Tau_k<-Tau_n*(!(model$q97.5$Tau>0 & model$q2.5$Tau<0))

par(mfrow=c(2,4),oma = c(3, 1, 2, 1))
cols = colorRampPalette(c("dark blue","white","red"))
col2 <- colorRampPalette(c("#4393C3", "#2166AC", "#053061",
                           "#FDDBC7", "#FFFFFF", "#D1E5F0", "#92C5DE",
                           "#67001F", "#B2182B", "#D6604D", "#F4A582"))


corrplot(Y_cor, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Correlation cor(Y)")
corrplot(me5$mean$EnvRho, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("EnvRho")
corrplot(me5$mean$EnvRho*(!me5$overlap0$EnvRho), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("EnvRho signif")
corrplot(me5$mean$Rho, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Rho")
corrplot(me5$mean$Rho*(!me5$overlap0$Rho), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Rho signif")
corrplot(Tau_n, diag = FALSE, order ="original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Tau")
corrplot(Tau_n*(!me5$overlap0$Tau), diag = FALSE, order ="original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Tau signif")

GJAM

##to setup chains parameters
data<-sim_data$EnvEvenSp5
#fit_gjam(data,2000,1000,"./gjam_models/gjam5env.rda")
load_gjam(data,name="./gjam_models/gjam5env.rda")
## 
## Sensitivity by predictor variables f:
##      Estimate   SE CI_025 CI_975
## env       9.5 3.13   5.08   16.4
## env2     31.9 7.93  20.50   50.6
## 
## Coefficient matrix B:
##              sp01   sp02    sp03   sp04   sp05
## intercept -0.0312  0.646  1.2500  0.704 -1.200
## env       -2.4200 -1.800  0.0438  1.710  2.870
## env2       0.9440 -0.527 -0.9410 -0.598 -0.665
## RMSPE      0.2780  0.327  0.3760  0.309  0.327
## 
## Coefficient matrix B:
##                Estimate     SE  CI_025 CI_975 sig95
## sp01_intercept  -0.0312 0.0775 -0.1940  0.107      
## sp01_env        -2.4200 0.1170 -2.6500 -2.200     *
## sp01_env2        0.9440 0.0938  0.7780  1.130     *
## sp02_intercept   0.6460 0.0651  0.5220  0.771     *
## sp02_env        -1.8000 0.0911 -1.9900 -1.630     *
## sp02_env2       -0.5270 0.0741 -0.6700 -0.389     *
## sp03_intercept   1.2500 0.0672  1.1100  1.380     *
## sp03_env         0.0438 0.0465 -0.0447  0.132      
## sp03_env2       -0.9410 0.0764 -1.0900 -0.781     *
## sp04_intercept   0.7040 0.0865  0.5590  0.881     *
## sp04_env         1.7100 0.0897  1.5500  1.890     *
## sp04_env2       -0.5980 0.0608 -0.7170 -0.486     *
## sp05_intercept  -1.2000 0.0953 -1.3700 -1.000     *
## sp05_env         2.8700 0.1380  2.6300  3.170     *
## sp05_env2       -0.6650 0.0610 -0.7840 -0.545     *
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X and W:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
##  Design Table
##      env env2
## VIF    1    1
## env2   0   NA
## 
##  Sample contains n = 500 observations on S = 5 response variables, and 2 predictors.  Data types (typeNames) include PA. There are 0 missing values in X and 0 missing values in Y. The RMSPE is 0.325, and the DIC is 25294.  Computation involved 2000 Gibbs steps, with a burnin of 1000.

HMSC

#setwd("~/Tesi/Code/Ecology-models-master/simcoms-master/ExampleFiles")
fit_hmsc<-function(data,label="Fit",nsamples = 1000,nchains=2,name="./HMmodels/hmtemp.rda" ){
  if (label=="Fit"){
    Y_data = subset(data, select = -env)
    ns<- ncol(Y_data)
    np <- nrow(Y_data)
    X<-scale(poly(data$env[1:np], 2))
    colnames(X)<-c("env","env2")
    studyDesign = data.frame(sample = as.factor(1:np))
    rL = HmscRandomLevel(units = studyDesign$sample)
    m = Hmsc(Y=as.matrix(Y_data), XData=as.data.frame(X), XFormula=~env+env2, distr="probit",
             studyDesign = studyDesign, ranLevels = list(sample = rL))
    m = sampleMcmc(m, nsamples, thin=10, adaptNf=c(200,200), transient=500,nChains=nchains ,verbose=F)
    save(m, file = name)
    return(m)
  }
  if (label=="Load"){
    return(load_object(name))
  }
}



data<-sim_data$EnvEvenSp5
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm5env.rda" )
#hm_mod<-load_object("./HMmodels/hm5env.rda")

Convergence:

hm_conv<-function(mod){
  codaList = convertToCodaObject(mod)

  #convergence histograms
  hist(effectiveSize(codaList$Beta), main="ess(beta)",lwd=2,col=gray(.6))
  hist(gelman.diag(codaList$Beta,multivariate=FALSE)$psrf,lwd=2,col=gray(.6), main="psrf(beta)")
  
  hist(effectiveSize(codaList$Omega[[1]]), main="ess(omega)",lwd=2,col=gray(.6))
  hist(gelman.diag(codaList$Omega[[1]], multivariate=FALSE)$psrf,lwd=2,col=gray(.6), main="psrf(omega)")
}

hm_conv(hm_mod)

Study of interactions

hm_inter<-function(mod, nchains=2,nsamples = 1000, interact=diag(ns)){
  getOmega = function(a,r=1)
  return(crossprod(a$Lambda[[r]]))
  ns<-mod$ns
  postOmega1 = array(unlist(lapply(mod$postList[[1]],getOmega)),c(ns,ns,mod$samples))
  postOmega2 = array(unlist(lapply(mod$postList[[2]],getOmega)),c(ns,ns,mod$samples))
  
  postOmega<-abind(postOmega1,postOmega2,along=3)
  postOmegaMean = apply(postOmega,c(1,2),mean)
  postOmegaUp=apply(postOmega,c(1,2),quantile,0.95)
  postOmegaLo=apply(postOmega,c(1,2),quantile,0.05)
  
  postR<-array(dim=c(ns,ns,nchains*nsamples))
  for(i in 1:dim(postOmega)[3])
  postR[,,i]<-stats::cov2cor(postOmega[,,i])
  
  
  postRMean = apply(postR,c(1,2),mean)
  postRUp=apply(postR,c(1,2),quantile,0.95)
  postRLo=apply(postR,c(1,2),quantile,0.05)
  
  Tau = solve(postOmegaMean)
  Tau_n = -cov2cor(Tau)
  
  
  Toplot_R<-postRMean*(!(postRUp>0 & postRLo<0))
  
  # Omegacor<- computeAssociations(m)
  # supportLevel<- 0.95
  # toPlot<- ((Omegacor[[1]]$support>supportLevel)+ (Omegacor[[1]]$support<(1-supportLevel))>0)*Omegacor[[1]]$mean
  # corrplot(toPlot, method="color", col=colorRampPalette(c("blue", "white", "red"))(200))
  
  par(mfrow=c(2,3),oma = c(1, 1, 1, 1))
  corrplot(cor(hm_mod$Y), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("Correlation cor(Y)")
  corrplot(postRMean, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("R")
  corrplot(Toplot_R, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("Plot only non zero value")
  corrplot(Tau_n, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("Partial correlation matrix")
  corrplot(interact, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("True interactions")
}

hm_inter(hm_mod)

Environmental filtering 10 species

JSDM

me10 <- load_object("model-2019-04-10-08-26-20.rda")
#load("model-2019-04-09-19-02-16.rda")
summary(me10)
## Summary for model '/tmp/RtmpKix4lZ/file40e957831178'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 683.734 minutes at time 2019-04-09 21:02:35.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
jsdm_conv(me10)
## Maximum Rhat value for Beta: 1.5602
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.7953

me10$mcmc.info[1:7]
## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
data<-sim_data$EnvEvenSp10
data <- list(
  Y = subset(data, select = -env),
  X = cbind(1, scale(poly(data$env, 2))),
  covx = cov(cbind(1, scale(poly(data$env, 2)))),
  K = 3,
  J = ncol(data) - 1,
  n = nrow(data),
  I = diag(ncol(data) - 1),
  df = ncol(data)
)

##########################################################################################
#Tau<-solve
#Tau_n<-to_prec(me10$mean$Tau)
#Tau_k<-Tau_n*(!(model$q97.5$Tau>0 & model$q2.5$Tau<0))

plot_cor_jsdm<-function(mod,y,interact=diag(ncol(y))){
  par(mfrow=c(2,4),oma = c(3, 1, 2, 1))
  corrplot(cor(y), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("Correlation cor(Y)")
  corrplot(mod$mean$EnvRho, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("EnvRho")
  corrplot(mod$mean$EnvRho*(!mod$overlap0$EnvRho), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("EnvRho signif")
  corrplot(mod$mean$Rho, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("Rho")
  corrplot(mod$mean$Rho*(!mod$overlap0$Rho), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("Rho signif")
  corrplot(interact, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
  title("True interactions")
}

plot_cor_jsdm(me10,data$Y)

Gjam

data<-sim_data$EnvEvenSp10

#fit_gjam(data,5000,500,"./gjam_models/gjam10env.rda")
load_gjam(data,name="./gjam_models/gjam10env.rda")
## 
## Sensitivity by predictor variables f:
##      Estimate   SE CI_025 CI_975
## env     233.0 58.0    120  354.0
## env2     56.7 11.8     37   82.9
## 
## Coefficient matrix B:
##              sp01    sp02   sp03   sp04   sp05   sp06   sp07   sp08  sp09
## intercept -0.2060  0.3800  1.310  0.889  1.020  0.742  0.638  0.140 0.332
## env       -2.0500 -2.5700 -2.300 -1.030 -0.189  0.542  1.090  1.740 2.520
## env2       0.0656 -0.0531  0.338 -0.533 -0.805 -0.838 -0.690 -0.425 0.640
## RMSPE      0.3300  0.2980  0.317  0.371  0.414  0.421  0.374  0.345 0.298
##             sp10
## intercept -0.577
## env        2.370
## env2       0.677
## RMSPE      0.276
## 
## Coefficient matrix B:
##                Estimate     SE   CI_025  CI_975 sig95
## sp01_intercept  -0.2060 0.0670 -0.33700 -0.0733     *
## sp01_env        -2.0500 0.1190 -2.29000 -1.8300     *
## sp01_env2        0.0656 0.0992 -0.14100  0.2380      
## sp02_intercept   0.3800 0.1030  0.18000  0.5510     *
## sp02_env        -2.5700 0.1370 -2.83000 -2.2900     *
## sp02_env2       -0.0531 0.0790 -0.19700  0.1120      
## sp03_intercept   1.3100 0.0829  1.15000  1.4700     *
## sp03_env        -2.3000 0.1260 -2.55000 -2.0600     *
## sp03_env2        0.3380 0.0712  0.20500  0.4860     *
## sp04_intercept   0.8890 0.0585  0.77400  1.0000     *
## sp04_env        -1.0300 0.0659 -1.16000 -0.9070     *
## sp04_env2       -0.5330 0.0580 -0.64200 -0.4160     *
## sp05_intercept   1.0200 0.0621  0.90800  1.1500     *
## sp05_env        -0.1890 0.0488 -0.28500 -0.0949     *
## sp05_env2       -0.8050 0.0759 -0.96000 -0.6590     *
## sp06_intercept   0.7420 0.0635  0.61600  0.8650     *
## sp06_env         0.5420 0.0558  0.43200  0.6460     *
## sp06_env2       -0.8380 0.0634 -0.96100 -0.7130     *
## sp07_intercept   0.6380 0.0795  0.50600  0.8160     *
## sp07_env         1.0900 0.0805  0.93900  1.2500     *
## sp07_env2       -0.6900 0.0595 -0.80300 -0.5700     *
## sp08_intercept   0.1400 0.0759  0.00798  0.2950     *
## sp08_env         1.7400 0.1240  1.48000  1.9600     *
## sp08_env2       -0.4250 0.1000 -0.65400 -0.2750     *
## sp09_intercept   0.3320 0.0539  0.22600  0.4370     *
## sp09_env         2.5200 0.1660  2.22000  2.8500     *
## sp09_env2        0.6400 0.0700  0.50400  0.7780     *
## sp10_intercept  -0.5770 0.0696 -0.70200 -0.4380     *
## sp10_env         2.3700 0.1130  2.16000  2.6000     *
## sp10_env2        0.6770 0.0820  0.49900  0.8200     *
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X and W:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
##  Design Table
##      env env2
## VIF    1    1
## env2   0   NA
## 
##  Sample contains n = 500 observations on S = 10 response variables, and 2 predictors.  Data types (typeNames) include PA. There are 0 missing values in X and 0 missing values in Y. The RMSPE is 0.348, and the DIC is 114934.  Computation involved 5000 Gibbs steps, with a burnin of 500.

#gje10<-load_object("./gjam_models/gjam10env.rda")

#to check posterior density of s in Sigma 
#gje10<-load_object("./gjam_models/gjam10env.rda")
#plot(density(gje10$chains$sgibbs[,4]))

HMSC

data<-sim_data$EnvEvenSp10
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm10env.rda" )

hm_conv(hm_mod)

hm_inter(hm_mod)

Environmental filtering 20 species

JSDM

me20 <- load_object("model-2019-04-11-19-06-02.rda")
#load("model-2019-04-09-19-02-16.rda")
summary(me20)
## Summary for model '/tmp/RtmpKix4lZ/file40e966e51ba7'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 2079.661 minutes at time 2019-04-10 08:26:21.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
jsdm_conv(me20)
## Maximum Rhat value for Beta: 1.2174
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.3236

me20$mcmc.info[1:7]
## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
data<-sim_data$EnvEvenSp20
data <- list(
  Y = subset(data, select = -env),
  X = cbind(1, scale(poly(data$env, 2))),
  covx = cov(cbind(1, scale(poly(data$env, 2)))),
  K = 3,
  J = ncol(data) - 1,
  n = nrow(data),
  I = diag(ncol(data) - 1),
  df = ncol(data)
)

##########################################################################################

plot_cor_jsdm(me20,data$Y)

Gjam

data<-sim_data$EnvEvenSp20

#fit_gjam(data,5000,500,"./gjam_models/gjam20env.rda")

load_gjam(data,name="./gjam_models/gjam20env.rda")
## 
## Sensitivity by predictor variables f:
##      Estimate   SE CI_025 CI_975
## env       335 61.7    222    461
## env2      157 28.3    110    219
## 
## Coefficient matrix B:
##             sp01    sp02   sp03   sp04   sp05   sp06   sp07   sp08   sp09
## intercept -0.533  0.0688  0.500  0.144  0.327  0.498  0.720  0.896  0.906
## env       -1.850 -2.5200 -2.480 -2.020 -1.820 -1.760 -1.220 -0.753 -0.516
## env2       0.187  0.3400  0.451 -0.279 -0.507 -0.493 -0.651 -0.641 -0.815
## RMSPE      0.337  0.3000  0.309  0.331  0.332  0.333  0.365  0.398  0.406
##             sp10   sp11   sp12   sp13   sp14   sp15   sp16   sp17     sp18
## intercept  1.130  1.120  1.050  0.689  0.610  0.435  0.440  0.113 -0.07100
## env       -0.200  0.133  0.551  0.870  1.140  1.730  2.460  1.990  2.35000
## env2      -0.931 -0.855 -0.860 -0.794 -0.671 -0.452 -0.227 -0.271  0.00271
## RMSPE      0.398  0.411  0.386  0.399  0.363  0.339  0.308  0.330  0.30300
##              sp19   sp20
## intercept -0.0697 -0.329
## env        2.5800  2.190
## env2       0.3120  0.530
## RMSPE      0.3030  0.318
## 
## Coefficient matrix B:
##                Estimate     SE    CI_025  CI_975 sig95
## sp01_intercept -0.53300 0.1340 -0.764000 -0.3250     *
## sp01_env       -1.85000 0.1340 -2.100000 -1.5800     *
## sp01_env2       0.18700 0.0717  0.055600  0.3320     *
## sp02_intercept  0.06880 0.0972 -0.095900  0.2760      
## sp02_env       -2.52000 0.1560 -2.810000 -2.2200     *
## sp02_env2       0.34000 0.0720  0.206000  0.4840     *
## sp03_intercept  0.50000 0.1040  0.295000  0.6900     *
## sp03_env       -2.48000 0.1360 -2.770000 -2.2400     *
## sp03_env2       0.45100 0.0813  0.310000  0.6170     *
## sp04_intercept  0.14400 0.0783  0.009940  0.3210     *
## sp04_env       -2.02000 0.0913 -2.200000 -1.8400     *
## sp04_env2      -0.27900 0.0825 -0.432000 -0.1160     *
## sp05_intercept  0.32700 0.0843  0.172000  0.4850     *
## sp05_env       -1.82000 0.1010 -2.020000 -1.6200     *
## sp05_env2      -0.50700 0.0716 -0.653000 -0.3720     *
## sp06_intercept  0.49800 0.0786  0.342000  0.6400     *
## sp06_env       -1.76000 0.0786 -1.920000 -1.6100     *
## sp06_env2      -0.49300 0.0669 -0.631000 -0.3680     *
## sp07_intercept  0.72000 0.0628  0.598000  0.8430     *
## sp07_env       -1.22000 0.0844 -1.370000 -1.0500     *
## sp07_env2      -0.65100 0.0816 -0.819000 -0.4980     *
## sp08_intercept  0.89600 0.0702  0.766000  1.0400     *
## sp08_env       -0.75300 0.0678 -0.886000 -0.6230     *
## sp08_env2      -0.64100 0.0614 -0.759000 -0.5180     *
## sp09_intercept  0.90600 0.0712  0.776000  1.0500     *
## sp09_env       -0.51600 0.0835 -0.663000 -0.3500     *
## sp09_env2      -0.81500 0.0669 -0.942000 -0.6790     *
## sp10_intercept  1.13000 0.0655  1.000000  1.2600     *
## sp10_env       -0.20000 0.0573 -0.305000 -0.0862     *
## sp10_env2      -0.93100 0.0676 -1.060000 -0.8010     *
## sp11_intercept  1.12000 0.0684  0.987000  1.2500     *
## sp11_env        0.13300 0.0565  0.022800  0.2410     *
## sp11_env2      -0.85500 0.0630 -0.978000 -0.7310     *
## sp12_intercept  1.05000 0.0730  0.906000  1.2000     *
## sp12_env        0.55100 0.0747  0.420000  0.7010     *
## sp12_env2      -0.86000 0.0752 -1.010000 -0.7180     *
## sp13_intercept  0.68900 0.0591  0.570000  0.8000     *
## sp13_env        0.87000 0.0813  0.709000  1.0300     *
## sp13_env2      -0.79400 0.0668 -0.922000 -0.6560     *
## sp14_intercept  0.61000 0.0641  0.481000  0.7300     *
## sp14_env        1.14000 0.0675  1.020000  1.2800     *
## sp14_env2      -0.67100 0.0665 -0.800000 -0.5470     *
## sp15_intercept  0.43500 0.0588  0.325000  0.5530     *
## sp15_env        1.73000 0.1250  1.520000  1.9800     *
## sp15_env2      -0.45200 0.0641 -0.579000 -0.3340     *
## sp16_intercept  0.44000 0.0521  0.337000  0.5420     *
## sp16_env        2.46000 0.1050  2.260000  2.6700     *
## sp16_env2      -0.22700 0.0501 -0.323000 -0.1280     *
## sp17_intercept  0.11300 0.0563  0.000618  0.2230     *
## sp17_env        1.99000 0.1310  1.760000  2.2600     *
## sp17_env2      -0.27100 0.0645 -0.394000 -0.1480     *
## sp18_intercept -0.07100 0.0549 -0.181000  0.0357      
## sp18_env        2.35000 0.1040  2.140000  2.5500     *
## sp18_env2       0.00271 0.0659 -0.130000  0.1260      
## sp19_intercept -0.06970 0.0875 -0.241000  0.0747      
## sp19_env        2.58000 0.1220  2.350000  2.8300     *
## sp19_env2       0.31200 0.0517  0.216000  0.4120     *
## sp20_intercept -0.32900 0.0556 -0.433000 -0.2160     *
## sp20_env        2.19000 0.0953  2.010000  2.3800     *
## sp20_env2       0.53000 0.0616  0.409000  0.6490     *
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X and W:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
##  Design Table
##      env env2
## VIF    1    1
## env2   0   NA
## 
##  Sample contains n = 500 observations on S = 20 response variables, and 2 predictors.  Data types (typeNames) include PA. There are 0 missing values in X and 0 missing values in Y. The RMSPE is 0.351, and the DIC is 371075.  Computation involved 5000 Gibbs steps, with a burnin of 500.

#gje20<-load_object("./gjam_models/gjam20env.rda")

#to check posterior density of s in Sigma 
#gje20<-load_object("./gjam_models/gjam20env.rda")
#plot(density(gje20$chains$sgibbs[,4]))

Gjam dimension reduction

data<-sim_data$EnvEvenSp20

#fit_gjam(data,5000,500,"./gjam_models/gjam20env.rda")

#load_gjam(data,name="./gjam_models/gjam20env.rda")
#gje20<-load_object("./gjam_models/gjam20env.rda")

#to check posterior density of s in Sigma 
#gje20<-load_object("./gjam_models/gjam20env.rda")
#plot(density(gje20$chains$sgibbs[,4]))
 data <- list(
    Y = subset(data, select = -env),
    X = cbind(1, scale(poly(data$env, 2))),
    covx = cov(cbind(1, scale(poly(data$env, 2)))),
    K = 3,
    J = ncol(data) - 1,
    n = nrow(data),
    I = diag(ncol(data) - 1),
    df = ncol(data)
  )
xdata<-as.data.frame(data$X[,-1])
colnames(xdata)<- c("env","env2")
ydata<-as.data.frame(data$Y)
###formula
rl   <- list(r = 8, N = 20)
formula<-as.formula( ~env+ env2)
ml  <- list(ng = 2500, burnin = 500, typeNames = 'PA', reductList = rl)
####fit


mod_gjam1  <- gjam(formula, xdata = xdata, ydata = ydata, modelList = ml)
## 
## Observations and responses:
## [1] 500  20
## Warning in .setupReduct(modelList, S, Q, n): dimension reduction requires
## reductList$N < no. responses
## 
## Dimension reduced from 20 X 20 -> 20 X 8 responses
## ===========================================================================
## expanding covariance chains
## ===========================================================================
save(mod_gjam1, file = "./gjam_models/gjam20env_dr.rda")
summary(mod_gjam1)
## 
## Sensitivity by predictor variables f:
##      Estimate     SE CI_025 CI_975
## env      1.12 0.0522   1.02   1.23
## env2     1.31 0.0608   1.19   1.43
## 
## Coefficient matrix B:
##             sp01    sp02      sp03   sp04   sp05   sp06   sp07   sp08
## intercept -0.331 -0.0900  0.023700  0.126  0.305  0.352  0.472  0.481
## env       -0.476 -0.4830 -0.482000 -0.477 -0.472 -0.465 -0.375 -0.433
## env2       0.280  0.0941 -0.000526 -0.110 -0.264 -0.327 -0.436 -0.456
## RMSPE      0.405  0.4040  0.407000  0.410  0.402  0.402  0.412  0.419
##             sp09   sp10   sp11   sp12   sp13   sp14   sp15   sp16   sp17
## intercept  0.485  0.489  0.489  0.487  0.475  0.478  0.378  0.252  0.153
## env       -0.430 -0.106  0.321  0.463  0.448  0.401  0.469  0.479  0.480
## env2      -0.478 -0.484 -0.475 -0.481 -0.478 -0.455 -0.346 -0.256 -0.167
## RMSPE      0.429  0.441  0.449  0.428  0.423  0.410  0.406  0.406  0.408
##              sp18   sp19   sp20
## intercept -0.0196 -0.117 -0.301
## env        0.4830  0.481  0.475
## env2       0.0111  0.104  0.246
## RMSPE      0.4030  0.406  0.403
## 
## Coefficient matrix B:
##                 Estimate     SE  CI_025  CI_975 sig95
## sp01_intercept -0.331000 0.0699 -0.4640 -0.1830     *
## sp01_env       -0.476000 0.0253 -0.4990 -0.4070     *
## sp01_env2       0.280000 0.0683  0.1400  0.4080     *
## sp02_intercept -0.090000 0.0703 -0.2290  0.0494      
## sp02_env       -0.483000 0.0170 -0.5000 -0.4400     *
## sp02_env2       0.094100 0.0702 -0.0424  0.2330      
## sp03_intercept  0.023700 0.0720 -0.1160  0.1680      
## sp03_env       -0.482000 0.0188 -0.5000 -0.4280     *
## sp03_env2      -0.000526 0.0686 -0.1320  0.1310      
## sp04_intercept  0.126000 0.0679 -0.0115  0.2540      
## sp04_env       -0.477000 0.0233 -0.4990 -0.4110     *
## sp04_env2      -0.110000 0.0690 -0.2490  0.0219      
## sp05_intercept  0.305000 0.0714  0.1670  0.4470     *
## sp05_env       -0.472000 0.0273 -0.4990 -0.4010     *
## sp05_env2      -0.264000 0.0721 -0.4070 -0.1240     *
## sp06_intercept  0.352000 0.0666  0.2200  0.4730     *
## sp06_env       -0.465000 0.0328 -0.4990 -0.3790     *
## sp06_env2      -0.327000 0.0669 -0.4550 -0.1900     *
## sp07_intercept  0.472000 0.0266  0.4000  0.4990     *
## sp07_env       -0.375000 0.0858 -0.4980 -0.2080     *
## sp07_env2      -0.436000 0.0471 -0.4970 -0.3270     *
## sp08_intercept  0.481000 0.0185  0.4290  0.4990     *
## sp08_env       -0.433000 0.0536 -0.4980 -0.2990     *
## sp08_env2      -0.456000 0.0364 -0.4990 -0.3660     *
## sp09_intercept  0.485000 0.0151  0.4440  0.5000     *
## sp09_env       -0.430000 0.0539 -0.4970 -0.3020     *
## sp09_env2      -0.478000 0.0206 -0.4990 -0.4230     *
## sp10_intercept  0.489000 0.0115  0.4560  0.5000     *
## sp10_env       -0.106000 0.0861 -0.2490  0.1060      
## sp10_env2      -0.484000 0.0157 -0.5000 -0.4420     *
## sp11_intercept  0.489000 0.0114  0.4600  0.5000     *
## sp11_env        0.321000 0.0707  0.1720  0.4540     *
## sp11_env2      -0.475000 0.0220 -0.4990 -0.4160     *
## sp12_intercept  0.487000 0.0129  0.4540  0.5000     *
## sp12_env        0.463000 0.0341  0.3710  0.4990     *
## sp12_env2      -0.481000 0.0175 -0.4990 -0.4350     *
## sp13_intercept  0.475000 0.0236  0.4100  0.4990     *
## sp13_env        0.448000 0.0495  0.3200  0.4990     *
## sp13_env2      -0.478000 0.0216 -0.5000 -0.4200     *
## sp14_intercept  0.478000 0.0208  0.4220  0.4990     *
## sp14_env        0.401000 0.0619  0.2730  0.4950     *
## sp14_env2      -0.455000 0.0383 -0.4990 -0.3570     *
## sp15_intercept  0.378000 0.0661  0.2360  0.4880     *
## sp15_env        0.469000 0.0279  0.3980  0.4990     *
## sp15_env2      -0.346000 0.0674 -0.4720 -0.2150     *
## sp16_intercept  0.252000 0.0670  0.1220  0.3810     *
## sp16_env        0.479000 0.0202  0.4230  0.4990     *
## sp16_env2      -0.256000 0.0679 -0.3850 -0.1240     *
## sp17_intercept  0.153000 0.0681  0.0191  0.2820     *
## sp17_env        0.480000 0.0205  0.4230  0.4990     *
## sp17_env2      -0.167000 0.0678 -0.2990 -0.0378     *
## sp18_intercept -0.019600 0.0713 -0.1580  0.1240      
## sp18_env        0.483000 0.0172  0.4370  0.5000     *
## sp18_env2       0.011100 0.0713 -0.1280  0.1490      
## sp19_intercept -0.117000 0.0701 -0.2560  0.0199      
## sp19_env        0.481000 0.0187  0.4310  0.5000     *
## sp19_env2       0.104000 0.0683 -0.0272  0.2420      
## sp20_intercept -0.301000 0.0693 -0.4320 -0.1640     *
## sp20_env        0.475000 0.0239  0.4120  0.4990     *
## sp20_env2       0.246000 0.0702  0.1100  0.3820     *
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X and W:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
##  Design Table
##      env env2
## VIF    1    1
## env2   0   NA
## 
##  Sample contains n = 500 observations on S = 20 response variables, and 2 predictors.  Data types (typeNames) include PA. There are 0 missing values in X and 0 missing values in Y. The RMSPE is 0.414, and the DIC is 77425.  Computation involved 2500 Gibbs steps, with a burnin of 500.  Dimension reduction was implemented with N = 9 and r = 8.
Tau <- solve(mod_gjam1$parameters$sigMu)
Tau_n = to_prec(Tau)

#postH<-apply(mod_gjam1$chains$sgibbs, 2, quantile,0.95)
#postL<-apply(mod_gjam1$chains$sgibbs, 2, quantile,0.05)


#pH<-convert_to_m(postH)
#pL<-convert_to_m(postL)

#R_sign<-cov2cor(mod_gjam1$parameters$sigMu)*(!(pH>0 & pL<0))

par(mfrow=c(2,3),oma = c(1, 1, 1, 1))
corrplot(cor(ydata), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Correlation cor(Y)")
corrplot(mod_gjam1$parameters$corMu, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("R")
corrplot(mod_gjam1$parameters$ematrix, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("E matrix")
# corrplot(Tau_n, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
# title("Tau")
# corrplot(R_sign, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
# title("R signif")
corrplot(diag(20), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("True interactions")

data<-sim_data$EnvEvenSp20
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm20env.rda" )
hm_conv(hm_mod)

hm_inter(hm_mod)

Environmental filtering + Facilitation dense 5 species

JSDM

mf5 <- load_object("model-2019-04-11-19-35-11.rda")
summary(mf5)
## Summary for model '/tmp/RtmpKix4lZ/file40e94e482135'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 29.13 minutes at time 2019-04-11 19:06:03.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
#mf5$Rhat
jsdm_conv(mf5)
## Maximum Rhat value for Beta: 1.3081
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.1421

mf5$mcmc.info[1:7]
## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
data<-sim_data$FacDenseSp5
data <- list(
  Y = subset(data, select = -env),
  X = cbind(1, scale(poly(data$env, 2))),
  covx = cov(cbind(1, scale(poly(data$env, 2)))),
  K = 3,
  J = ncol(data) - 1,
  n = nrow(data),
  I = diag(ncol(data) - 1),
  df = ncol(data)
)

##########################################################################################
#Tau<-solve
#Tau_n<-to_prec(me10$mean$Tau)
#Tau_k<-Tau_n*(!(model$q97.5$Tau>0 & model$q2.5$Tau<0))

plot_cor_jsdm(mf5,data$Y,fac_inter[[4]])

Gjam

data<-sim_data$FacDenseSp5
#fit_gjam(data,10000,2000,"./gjam_models/gjam5f.rda",interact=fac_inter[[4]])
load_gjam(data,name="./gjam_models/gjam5f.rda", interact=fac_inter[[4]])
## 
## Sensitivity by predictor variables f:
##      Estimate    SE CI_025 CI_975
## env      27.6 13.50   8.15   58.8
## env2     18.4  6.44   8.83   34.3
## 
## Coefficient matrix B:
##             sp01   sp02   sp03  sp04   sp05
## intercept -1.640  0.979  3.170 1.890 -1.840
## env       -1.750 -1.580  0.253 2.970  1.740
## env2       0.295 -0.575 -1.280 0.446  0.289
## RMSPE      0.310  0.312  0.204 0.326  0.308
## 
## Coefficient matrix B:
##                Estimate     SE CI_025 CI_975 sig95
## sp01_intercept   -1.640 0.1120 -1.910 -1.450     *
## sp01_env         -1.750 0.1450 -2.040 -1.490     *
## sp01_env2         0.295 0.0642  0.171  0.418     *
## sp02_intercept    0.979 0.1040  0.774  1.190     *
## sp02_env         -1.580 0.1390 -1.850 -1.330     *
## sp02_env2        -0.575 0.1090 -0.761 -0.347     *
## sp03_intercept    3.170 0.2530  2.730  3.680     *
## sp03_env          0.253 0.0773  0.110  0.402     *
## sp03_env2        -1.280 0.1260 -1.560 -1.050     *
## sp04_intercept    1.890 0.1980  1.450  2.230     *
## sp04_env          2.970 0.2390  2.500  3.450     *
## sp04_env2         0.446 0.1150  0.220  0.645     *
## sp05_intercept   -1.840 0.1670 -2.210 -1.570     *
## sp05_env          1.740 0.1260  1.500  1.990     *
## sp05_env2         0.289 0.0699  0.141  0.425     *
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X and W:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
##  Design Table
##      env env2
## VIF    1    1
## env2   0   NA
## 
##  Sample contains n = 500 observations on S = 5 response variables, and 2 predictors.  Data types (typeNames) include PA. There are 0 missing values in X and 0 missing values in Y. The RMSPE is 0.296, and the DIC is 45158.  Computation involved 10000 Gibbs steps, with a burnin of 2000.

#gjfd5<-load_object("./gjam_models/gjam5f.rda")

#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam5f.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))

HMSC

data<-sim_data$FacDenseSp5
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm5fd.rda" )
hm_conv(hm_mod)

hm_inter(hm_mod, interact = fac_inter[[4]])

Environmental filtering + Facilitation dense 10 species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e96d09e31e'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 684.487 minutes at time 2019-04-11 19:35:12.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.2866
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.4808

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000

data<-sim_data$FacDenseSp10

#fit_gjam(data,5000,500,"./gjam_models/gjam10fd.rda",interact=fac_inter[[5]])
load_gjam(data,name="./gjam_models/gjam10fd.rda", interact=fac_inter[[5]])
## 
## Sensitivity by predictor variables f:
##      Estimate    SE CI_025 CI_975
## env     160.0 35.60  101.0    240
## env2     18.5  4.88   10.4     29
## 
## Coefficient matrix B:
##             sp01   sp02   sp03   sp04    sp05   sp06   sp07   sp08  sp09
## intercept  0.212 -1.490 -0.877  2.410 -0.0951  2.400  0.497 -0.401 1.260
## env       -2.760 -2.660 -1.810 -3.010 -0.5910  1.160  1.330  1.030 3.990
## env2       1.070 -0.969 -0.644  0.319 -0.7020 -0.454 -0.256 -0.347 0.095
## RMSPE      0.274  0.375  0.406  0.305  0.4930  0.259  0.390  0.471 0.248
##             sp10
## intercept -0.893
## env        2.250
## env2      -0.411
## RMSPE      0.354
## 
## Coefficient matrix B:
##                Estimate     SE  CI_025  CI_975 sig95
## sp01_intercept   0.2120 0.0687  0.0833  0.3520     *
## sp01_env        -2.7600 0.1740 -3.1200 -2.4700     *
## sp01_env2        1.0700 0.1350  0.8190  1.3200     *
## sp02_intercept  -1.4900 0.2000 -1.9200 -1.1600     *
## sp02_env        -2.6600 0.2110 -3.0900 -2.3000     *
## sp02_env2       -0.9690 0.1130 -1.2400 -0.7780     *
## sp03_intercept  -0.8770 0.1070 -1.0900 -0.6900     *
## sp03_env        -1.8100 0.1470 -2.1200 -1.5700     *
## sp03_env2       -0.6440 0.1190 -0.9250 -0.4570     *
## sp04_intercept   2.4100 0.2040  2.0500  2.8000     *
## sp04_env        -3.0100 0.2470 -3.4600 -2.5900     *
## sp04_env2        0.3190 0.0696  0.1860  0.4540     *
## sp05_intercept  -0.0951 0.0554 -0.2050  0.0147      
## sp05_env        -0.5910 0.0628 -0.7140 -0.4700     *
## sp05_env2       -0.7020 0.0743 -0.8490 -0.5590     *
## sp06_intercept   2.4000 0.1260  2.1600  2.6400     *
## sp06_env         1.1600 0.0761  1.0100  1.3100     *
## sp06_env2       -0.4540 0.0649 -0.5840 -0.3310     *
## sp07_intercept   0.4970 0.0563  0.3840  0.6040     *
## sp07_env         1.3300 0.0797  1.1800  1.4900     *
## sp07_env2       -0.2560 0.0554 -0.3670 -0.1470     *
## sp08_intercept  -0.4010 0.0539 -0.5070 -0.2950     *
## sp08_env         1.0300 0.0877  0.8540  1.1900     *
## sp08_env2       -0.3470 0.0592 -0.4640 -0.2320     *
## sp09_intercept   1.2600 0.0922  1.0900  1.4500     *
## sp09_env         3.9900 0.3810  3.3800  4.6800     *
## sp09_env2        0.0950 0.0613 -0.0361  0.2050      
## sp10_intercept  -0.8930 0.0796 -1.0700 -0.7550     *
## sp10_env         2.2500 0.1280  2.0400  2.5200     *
## sp10_env2       -0.4110 0.0576 -0.5230 -0.2980     *
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X and W:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
##  Design Table
##      env env2
## VIF    1    1
## env2   0   NA
## 
##  Sample contains n = 500 observations on S = 10 response variables, and 2 predictors.  Data types (typeNames) include PA. There are 0 missing values in X and 0 missing values in Y. The RMSPE is 0.367, and the DIC is 157914.  Computation involved 5000 Gibbs steps, with a burnin of 500.

#gjfd5<-load_object("./gjam_models/gjam10fd.rda")

#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam10fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))

HMSC

Environmental filtering + Facilitation dense 20 species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e91ec9ff9'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 2062.467 minutes at time 2019-04-12 06:59:42.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.6132
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.508

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000

Gjam

data<-sim_data$FacDenseSp20

#fit_gjam(data,10000,1500,"./gjam_models/gjam20fd.rda",interact=fac_inter[[6]])
load_gjam(data,name="./gjam_models/gjam20fd.rda", interact=fac_inter[[6]])
## 
## Sensitivity by predictor variables f:
##      Estimate   SE CI_025 CI_975
## env       301 55.7    207    421
## env2      170 30.7    121    243
## 
## Coefficient matrix B:
##             sp01   sp02   sp03   sp04   sp05   sp06    sp07   sp08   sp09
## intercept -0.528 -1.010 -0.866  0.194 -0.280  0.327  0.0883  0.310  0.494
## env       -1.940 -1.950 -2.150 -2.820 -1.560 -2.170 -0.9060 -0.695 -0.778
## env2       0.294 -0.248 -0.681 -0.305 -0.784 -0.841 -0.6030 -0.903 -1.080
## RMSPE      0.323  0.362  0.371  0.293  0.401  0.285  0.4550  0.443  0.397
##             sp10   sp11   sp12   sp13   sp14   sp15   sp16  sp17    sp18
## intercept  0.696  0.396  0.315  0.525  0.118 -0.218 -0.186 0.336 0.00648
## env       -0.303  0.104  0.464  0.896  1.160  1.780  2.140 2.240 2.44000
## env2      -1.050 -0.883 -0.877 -0.915 -1.010 -1.260 -0.574 0.152 0.08790
## RMSPE      0.418  0.468  0.455  0.397  0.404  0.365  0.329 0.323 0.29800
##              sp19   sp20
## intercept -0.0434 -0.393
## env        3.1700  1.980
## env2       0.5820  0.270
## RMSPE      0.2720  0.320
## 
## Coefficient matrix B:
##                Estimate     SE   CI_025  CI_975 sig95
## sp01_intercept -0.52800 0.1010 -0.71400 -0.3420     *
## sp01_env       -1.94000 0.1660 -2.34000 -1.6800     *
## sp01_env2       0.29400 0.1980  0.02630  0.6830     *
## sp02_intercept -1.01000 0.1210 -1.25000 -0.7800     *
## sp02_env       -1.95000 0.1570 -2.25000 -1.6600     *
## sp02_env2      -0.24800 0.1100 -0.44500 -0.0388     *
## sp03_intercept -0.86600 0.2110 -1.22000 -0.4330     *
## sp03_env       -2.15000 0.1640 -2.44000 -1.7900     *
## sp03_env2      -0.68100 0.1130 -0.92200 -0.4750     *
## sp04_intercept  0.19400 0.0898  0.03490  0.3640     *
## sp04_env       -2.82000 0.1770 -3.17000 -2.4900     *
## sp04_env2      -0.30500 0.0689 -0.43400 -0.1700     *
## sp05_intercept -0.28000 0.1190 -0.48000 -0.0615     *
## sp05_env       -1.56000 0.1350 -1.80000 -1.3100     *
## sp05_env2      -0.78400 0.1100 -0.97000 -0.5580     *
## sp06_intercept  0.32700 0.0954  0.15300  0.5120     *
## sp06_env       -2.17000 0.1190 -2.41000 -1.9500     *
## sp06_env2      -0.84100 0.1240 -1.15000 -0.6270     *
## sp07_intercept  0.08830 0.0641 -0.03460  0.2160      
## sp07_env       -0.90600 0.0843 -1.07000 -0.7410     *
## sp07_env2      -0.60300 0.0704 -0.74300 -0.4710     *
## sp08_intercept  0.31000 0.0662  0.18200  0.4390     *
## sp08_env       -0.69500 0.0783 -0.84500 -0.5450     *
## sp08_env2      -0.90300 0.0732 -1.04000 -0.7580     *
## sp09_intercept  0.49400 0.0684  0.35900  0.6240     *
## sp09_env       -0.77800 0.0880 -0.94700 -0.6100     *
## sp09_env2      -1.08000 0.0683 -1.21000 -0.9440     *
## sp10_intercept  0.69600 0.0602  0.57300  0.8090     *
## sp10_env       -0.30300 0.0736 -0.44100 -0.1630     *
## sp10_env2      -1.05000 0.0715 -1.19000 -0.9050     *
## sp11_intercept  0.39600 0.0692  0.25800  0.5310     *
## sp11_env        0.10400 0.0612 -0.01520  0.2260      
## sp11_env2      -0.88300 0.0692 -1.01000 -0.7440     *
## sp12_intercept  0.31500 0.0579  0.20300  0.4270     *
## sp12_env        0.46400 0.0548  0.35400  0.5690     *
## sp12_env2      -0.87700 0.0644 -0.99800 -0.7500     *
## sp13_intercept  0.52500 0.0640  0.40400  0.6530     *
## sp13_env        0.89600 0.0644  0.76900  1.0200     *
## sp13_env2      -0.91500 0.0909 -1.08000 -0.7390     *
## sp14_intercept  0.11800 0.0485  0.02210  0.2120     *
## sp14_env        1.16000 0.0696  1.02000  1.2900     *
## sp14_env2      -1.01000 0.0635 -1.13000 -0.8850     *
## sp15_intercept -0.21800 0.0620 -0.33500 -0.0928     *
## sp15_env        1.78000 0.1120  1.57000  1.9900     *
## sp15_env2      -1.26000 0.0795 -1.42000 -1.1100     *
## sp16_intercept -0.18600 0.1770 -0.47300  0.1040      
## sp16_env        2.14000 0.1120  1.93000  2.3700     *
## sp16_env2      -0.57400 0.0868 -0.73500 -0.4050     *
## sp17_intercept  0.33600 0.0779  0.17800  0.4820     *
## sp17_env        2.24000 0.2570  1.87000  2.7400     *
## sp17_env2       0.15200 0.0798  0.00565  0.3220     *
## sp18_intercept  0.00648 0.0759 -0.13200  0.1440      
## sp18_env        2.44000 0.1220  2.20000  2.6900     *
## sp18_env2       0.08790 0.0509 -0.01750  0.1850      
## sp19_intercept -0.04340 0.0733 -0.16700  0.1070      
## sp19_env        3.17000 0.1380  2.91000  3.4400     *
## sp19_env2       0.58200 0.0714  0.44700  0.7220     *
## sp20_intercept -0.39300 0.0627 -0.50600 -0.2600     *
## sp20_env        1.98000 0.1130  1.78000  2.2200     *
## sp20_env2       0.27000 0.0674  0.14400  0.4060     *
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X and W:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
##  Design Table
##      env env2
## VIF    1    1
## env2   0   NA
## 
##  Sample contains n = 500 observations on S = 20 response variables, and 2 predictors.  Data types (typeNames) include PA. There are 0 missing values in X and 0 missing values in Y. The RMSPE is 0.374, and the DIC is 400202.  Computation involved 10000 Gibbs steps, with a burnin of 1500.

#gjfd5<-load_object("./gjam_models/gjam20fd.rda")

#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))

HMSC

data<-sim_data$FacDenseSp20
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm20fd.rda" )
hm_conv(hm_mod)

hm_inter(hm_mod, interact = fac_inter[[6]])

Environmental filtering + Facilitation sparse 5 species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e917d13dae'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 29.05 minutes at time 2019-04-13 17:22:13.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.373
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.4048

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000

Gjam

data<-sim_data$FacSparseSp5

#fit_gjam(data,2500,500,"./gjam_models/gjam5fs.rda",interact=fac_inter[[7]])
load_gjam(data,name="./gjam_models/gjam5fs.rda", interact=fac_inter[[7]])
## 
## Sensitivity by predictor variables f:
##      Estimate   SE CI_025 CI_975
## env      14.0 8.42   4.23   35.1
## env2      7.8 3.77   3.27   18.2
## 
## Coefficient matrix B:
##             sp01   sp02   sp03   sp04    sp05
## intercept -0.209  0.212  0.308  1.030  0.0554
## env       -2.850 -2.290 -0.671  2.060  4.3800
## env2       0.689 -0.895 -1.020 -0.513 -0.0713
## RMSPE      0.284  0.300  0.431  0.281  0.2490
## 
## Coefficient matrix B:
##                Estimate     SE  CI_025   CI_975 sig95
## sp01_intercept  -0.2090 0.0925 -0.3770 -0.00693     *
## sp01_env        -2.8500 0.2760 -3.4700 -2.45000     *
## sp01_env2        0.6890 0.1140  0.4510  0.90200     *
## sp02_intercept   0.2120 0.0582  0.0951  0.33100     *
## sp02_env        -2.2900 0.1260 -2.5300 -2.05000     *
## sp02_env2       -0.8950 0.0888 -1.0600 -0.72600     *
## sp03_intercept   0.3080 0.0629  0.1830  0.43100     *
## sp03_env        -0.6710 0.0731 -0.8110 -0.53600     *
## sp03_env2       -1.0200 0.0822 -1.1700 -0.85900     *
## sp04_intercept   1.0300 0.1100  0.8240  1.25000     *
## sp04_env         2.0600 0.1490  1.8200  2.39000     *
## sp04_env2       -0.5130 0.0616 -0.6320 -0.39400     *
## sp05_intercept   0.0554 0.1810 -0.2860  0.42500      
## sp05_env         4.3800 0.6540  3.3400  5.72000     *
## sp05_env2       -0.0713 0.1730 -0.3960  0.28900      
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X and W:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
##  Design Table
##      env env2
## VIF    1    1
## env2   0   NA
## 
##  Sample contains n = 500 observations on S = 5 response variables, and 2 predictors.  Data types (typeNames) include PA. There are 0 missing values in X and 0 missing values in Y. The RMSPE is 0.315, and the DIC is 18491.  Computation involved 2500 Gibbs steps, with a burnin of 500.

#gjfd5<-load_object("./gjam_models/gjam5fs.rda")

#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))

HMSC

data<-sim_data$FacSparseSp5

hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm5fs.rda" )
hm_conv(hm_mod)

hm_inter(hm_mod, interact = fac_inter[[7]])

Environmental filtering + Facilitation sparse 10 species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e9113ddc9'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 686.687 minutes at time 2019-04-13 17:51:16.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.5479
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.4746

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000

data<-sim_data$FacSparseSp10

#fit_gjam(data,2500,500,"./gjam_models/gjam10fs.rda",interact=fac_inter[[8]])
load_gjam(data,name="./gjam_models/gjam10fs.rda", interact=fac_inter[[8]])
## 
## Sensitivity by predictor variables f:
##      Estimate   SE CI_025 CI_975
## env     162.0 39.3   92.4  247.0
## env2     59.7 12.6   39.3   88.8
## 
## Coefficient matrix B:
##             sp01   sp02   sp03   sp04   sp05   sp06   sp07    sp08   sp09
## intercept -1.140  0.492  0.732  0.052  0.278  0.268  0.108  0.9400 -0.888
## env       -1.640 -2.630 -2.650 -0.992 -0.158  0.341  0.861  2.8900  1.840
## env2      -0.304  0.102 -0.413 -1.290 -1.130 -1.010 -0.984 -0.0437 -0.799
## RMSPE      0.412  0.297  0.267  0.408  0.434  0.451  0.436  0.2940  0.407
##              sp10
## intercept -0.5640
## env        3.5400
## env2      -0.0422
## RMSPE      0.2560
## 
## Coefficient matrix B:
##                Estimate     SE   CI_025  CI_975 sig95
## sp01_intercept  -1.1400 0.0911 -1.32000 -0.9660     *
## sp01_env        -1.6400 0.1090 -1.83000 -1.4100     *
## sp01_env2       -0.3040 0.0872 -0.44700 -0.1160     *
## sp02_intercept   0.4920 0.0847  0.32600  0.6620     *
## sp02_env        -2.6300 0.1210 -2.87000 -2.3900     *
## sp02_env2        0.1020 0.0618 -0.00800  0.2300      
## sp03_intercept   0.7320 0.0855  0.57700  0.8950     *
## sp03_env        -2.6500 0.1360 -2.94000 -2.4100     *
## sp03_env2       -0.4130 0.0611 -0.53900 -0.3000     *
## sp04_intercept   0.0520 0.0733 -0.08990  0.1830      
## sp04_env        -0.9920 0.0793 -1.17000 -0.8600     *
## sp04_env2       -1.2900 0.0841 -1.45000 -1.1200     *
## sp05_intercept   0.2780 0.0621  0.16600  0.4120     *
## sp05_env        -0.1580 0.0607 -0.27700 -0.0445     *
## sp05_env2       -1.1300 0.0726 -1.28000 -0.9980     *
## sp06_intercept   0.2680 0.0596  0.15700  0.3860     *
## sp06_env         0.3410 0.0632  0.21300  0.4600     *
## sp06_env2       -1.0100 0.0733 -1.15000 -0.8710     *
## sp07_intercept   0.1080 0.0503  0.00799  0.2060     *
## sp07_env         0.8610 0.0894  0.70400  1.0300     *
## sp07_env2       -0.9840 0.0698 -1.12000 -0.8490     *
## sp08_intercept   0.9400 0.0703  0.79400  1.0700     *
## sp08_env         2.8900 0.1680  2.54000  3.2100     *
## sp08_env2       -0.0437 0.0504 -0.13900  0.0564      
## sp09_intercept  -0.8880 0.0713 -1.02000 -0.7410     *
## sp09_env         1.8400 0.0811  1.67000  1.9900     *
## sp09_env2       -0.7990 0.0645 -0.91600 -0.6650     *
## sp10_intercept  -0.5640 0.0605 -0.67500 -0.4400     *
## sp10_env         3.5400 0.1800  3.23000  3.9300     *
## sp10_env2       -0.0422 0.0522 -0.14100  0.0606      
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X and W:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
##  Design Table
##      env env2
## VIF    1    1
## env2   0   NA
## 
##  Sample contains n = 500 observations on S = 10 response variables, and 2 predictors.  Data types (typeNames) include PA. There are 0 missing values in X and 0 missing values in Y. The RMSPE is 0.374, and the DIC is 120194.  Computation involved 2500 Gibbs steps, with a burnin of 500.

#gjfd5<-load_object("./gjam_models/gjam10fs.rda")

#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam10fs.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))

HMSC

data<-sim_data$FacSparseSp10

hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm10fs.rda" )
hm_conv(hm_mod)

hm_inter(hm_mod, interact = fac_inter[[8]])

Environmental filtering + Facilitation sparse 20 species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e9629c97ad'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 2056.316 minutes at time 2019-04-14 05:17:59.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.5653
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.4826

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000

Gjam

data<-sim_data$FacSparseSp20

#fit_gjam(data,10000,1500,"./gjam_models/gjam20fs.rda",interact=fac_inter[[9]])
load_gjam(data,name="./gjam_models/gjam20fs.rda", interact=fac_inter[[9]])
## 
## Sensitivity by predictor variables f:
##      Estimate   SE CI_025 CI_975
## env       300 65.9  188.0    443
## env2      148 31.5   94.4    218
## 
## Coefficient matrix B:
##             sp01   sp02   sp03    sp04    sp05   sp06   sp07   sp08   sp09
## intercept -0.330 -1.030  0.197 -0.0566  0.6180  0.691  0.485  0.555  0.880
## env       -2.590 -2.240 -2.560 -2.1000 -1.7500 -1.690 -1.130 -0.915 -0.641
## env2       0.429 -0.362  0.131 -0.5870 -0.0122 -0.439 -0.559 -0.728 -0.860
## RMSPE      0.295  0.348  0.300  0.3050  0.3520  0.333  0.395  0.410  0.403
##             sp10   sp11   sp12   sp13   sp14   sp15   sp16   sp17  sp18
## intercept  0.862  0.847  0.601  0.700  0.447  0.594  0.337  0.302 0.425
## env       -0.212  0.236  0.513  1.050  0.987  1.530  1.520  2.750 2.190
## env2      -0.988 -1.140 -0.824 -0.958 -0.741 -0.448 -0.281 -0.307 0.525
## RMSPE      0.410  0.401  0.438  0.363  0.408  0.341  0.364  0.277 0.320
##             sp19   sp20
## intercept -0.105 -0.724
## env        2.440  2.280
## env2       0.281 -0.001
## RMSPE      0.292  0.304
## 
## Coefficient matrix B:
##                Estimate     SE  CI_025  CI_975 sig95
## sp01_intercept  -0.3300 0.0646 -0.4540 -0.1970     *
## sp01_env        -2.5900 0.1480 -2.8600 -2.2700     *
## sp01_env2        0.4290 0.1770  0.0651  0.6740     *
## sp02_intercept  -1.0300 0.1010 -1.2200 -0.8380     *
## sp02_env        -2.2400 0.1410 -2.5600 -2.0100     *
## sp02_env2       -0.3620 0.1120 -0.6120 -0.1510     *
## sp03_intercept   0.1970 0.0816  0.0411  0.3590     *
## sp03_env        -2.5600 0.1810 -2.9200 -2.2200     *
## sp03_env2        0.1310 0.1050 -0.0655  0.3510      
## sp04_intercept  -0.0566 0.0891 -0.2220  0.1230      
## sp04_env        -2.1000 0.1340 -2.4000 -1.8700     *
## sp04_env2       -0.5870 0.0776 -0.7300 -0.4350     *
## sp05_intercept   0.6180 0.0920  0.4430  0.7790     *
## sp05_env        -1.7500 0.1070 -1.9600 -1.5400     *
## sp05_env2       -0.0122 0.1040 -0.2690  0.1470      
## sp06_intercept   0.6910 0.0806  0.5250  0.8370     *
## sp06_env        -1.6900 0.1130 -1.9200 -1.4900     *
## sp06_env2       -0.4390 0.1030 -0.6550 -0.2760     *
## sp07_intercept   0.4850 0.0750  0.3440  0.6410     *
## sp07_env        -1.1300 0.1010 -1.3300 -0.9460     *
## sp07_env2       -0.5590 0.1180 -0.7940 -0.3560     *
## sp08_intercept   0.5550 0.0615  0.4290  0.6730     *
## sp08_env        -0.9150 0.0719 -1.0600 -0.7820     *
## sp08_env2       -0.7280 0.0797 -0.8670 -0.5640     *
## sp09_intercept   0.8800 0.0866  0.6900  1.0300     *
## sp09_env        -0.6410 0.0609 -0.7610 -0.5230     *
## sp09_env2       -0.8600 0.1060 -1.0700 -0.6590     *
## sp10_intercept   0.8620 0.0572  0.7490  0.9730     *
## sp10_env        -0.2120 0.0564 -0.3190 -0.0966     *
## sp10_env2       -0.9880 0.0745 -1.1400 -0.8510     *
## sp11_intercept   0.8470 0.0636  0.7210  0.9710     *
## sp11_env         0.2360 0.0671  0.1140  0.3790     *
## sp11_env2       -1.1400 0.0851 -1.2900 -0.9680     *
## sp12_intercept   0.6010 0.0696  0.4510  0.7280     *
## sp12_env         0.5130 0.0584  0.3980  0.6270     *
## sp12_env2       -0.8240 0.0719 -0.9680 -0.6890     *
## sp13_intercept   0.7000 0.0739  0.5550  0.8400     *
## sp13_env         1.0500 0.0776  0.9050  1.2000     *
## sp13_env2       -0.9580 0.0740 -1.1000 -0.8100     *
## sp14_intercept   0.4470 0.0761  0.2960  0.5960     *
## sp14_env         0.9870 0.0819  0.8290  1.1500     *
## sp14_env2       -0.7410 0.0739 -0.8830 -0.5950     *
## sp15_intercept   0.5940 0.0578  0.4820  0.7060     *
## sp15_env         1.5300 0.0960  1.3500  1.7100     *
## sp15_env2       -0.4480 0.0678 -0.5700 -0.3110     *
## sp16_intercept   0.3370 0.0589  0.2160  0.4480     *
## sp16_env         1.5200 0.0753  1.3700  1.6700     *
## sp16_env2       -0.2810 0.0615 -0.4050 -0.1670     *
## sp17_intercept   0.3020 0.1020  0.1030  0.4720     *
## sp17_env         2.7500 0.1330  2.4700  2.9900     *
## sp17_env2       -0.3070 0.0872 -0.4540 -0.1390     *
## sp18_intercept   0.4250 0.0906  0.2660  0.6020     *
## sp18_env         2.1900 0.1370  1.9400  2.4600     *
## sp18_env2        0.5250 0.0868  0.3550  0.6790     *
## sp19_intercept  -0.1050 0.0792 -0.2810  0.0245      
## sp19_env         2.4400 0.1240  2.2200  2.7000     *
## sp19_env2        0.2810 0.0574  0.1650  0.3890     *
## sp20_intercept  -0.7240 0.0791 -0.8750 -0.5750     *
## sp20_env         2.2800 0.1140  2.0600  2.5000     *
## sp20_env2       -0.0010 0.0740 -0.1410  0.1360      
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X and W:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
##  Design Table
##      env env2
## VIF    1    1
## env2   0   NA
## 
##  Sample contains n = 500 observations on S = 20 response variables, and 2 predictors.  Data types (typeNames) include PA. There are 0 missing values in X and 0 missing values in Y. The RMSPE is 0.356, and the DIC is 450253.  Computation involved 10000 Gibbs steps, with a burnin of 1500.

#gjfd5<-load_object("./gjam_models/gjam20fd.rda")

#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))

HMSC

HMSC

data<-sim_data$FacSparseSp20

hm_mod<-fit_hmsc(data,nsamples=3000, nchains=2,"Load",name="./HMmodels/hm20fs.rda" )
hm_conv(hm_mod)

hm_inter(hm_mod, nsamples=3000, nchains=2,interact = fac_inter[[9]])

Environmental filtering + Competition dense 5 species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e9127bbb96'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 29.308 minutes at time 2019-04-15 15:34:20.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## Successful convergence based on Rhat values (all < 1.1).
## Maximum Rhat value for Beta: 1.0324
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.0313

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000

Gjam

data<-sim_data$CompDenseSp5

#fit_gjam(data,10000,1000,"./gjam_models/gjam5cmpd.rda",interact=comp_inter[[10]])

load_gjam(data,name="./gjam_models/gjam5cmpd.rda", interact=(-1)*comp_inter[[10]])
## 
## Sensitivity by predictor variables f:
##      Estimate   SE CI_025 CI_975
## env      73.8 28.4   22.3  132.0
## env2     43.7 19.0   17.4   89.3
## 
## Coefficient matrix B:
##             sp01   sp02    sp03   sp04   sp05
## intercept -1.340  0.170  1.7800  0.236 -1.420
## env       -0.661 -0.863  0.0249  0.662  0.680
## env2       0.280 -1.080 -0.7440 -0.949  0.287
## RMSPE      0.351  0.427  0.3080  0.456  0.343
## 
## Coefficient matrix B:
##                Estimate     SE  CI_025 CI_975 sig95
## sp01_intercept  -1.3400 0.0888 -1.5100 -1.170     *
## sp01_env        -0.6610 0.0777 -0.8020 -0.500     *
## sp01_env2        0.2800 0.0782  0.1170  0.417     *
## sp02_intercept   0.1700 0.0748  0.0112  0.302     *
## sp02_env        -0.8630 0.0689 -0.9970 -0.731     *
## sp02_env2       -1.0800 0.0847 -1.2500 -0.918     *
## sp03_intercept   1.7800 0.1300  1.5300  2.020     *
## sp03_env         0.0249 0.0622 -0.0989  0.147      
## sp03_env2       -0.7440 0.1000 -0.9310 -0.551     *
## sp04_intercept   0.2360 0.0543  0.1300  0.342     *
## sp04_env         0.6620 0.0765  0.5340  0.851     *
## sp04_env2       -0.9490 0.0790 -1.1200 -0.807     *
## sp05_intercept  -1.4200 0.0901 -1.6000 -1.240     *
## sp05_env         0.6800 0.0742  0.5440  0.831     *
## sp05_env2        0.2870 0.0746  0.1490  0.425     *
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X and W:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
##  Design Table
##      env env2
## VIF    1    1
## env2   0   NA
## 
##  Sample contains n = 500 observations on S = 5 response variables, and 2 predictors.  Data types (typeNames) include PA. There are 0 missing values in X and 0 missing values in Y. The RMSPE is 0.381, and the DIC is 31214.  Computation involved 10000 Gibbs steps, with a burnin of 1000.

#gjfd5<-load_object("./gjam_models/gjam5cmpd.rda")

#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam5cmpd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
data<-sim_data$CompDenseSp5

hm_mod<-fit_hmsc(data,"Load",nsamples=3000, nchains=2,name="./HMmodels/hm5cmpd.rda" )
hm_conv(hm_mod)

hm_inter(hm_mod, nsamples=3000, nchains=2,interact = (-1)*comp_inter[[10]])

Environmental filtering + Competition dense 10 species

## Summary for model '/tmp/RtmpKix4lZ/file40e9559ce9d8'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 695.614 minutes at time 2019-04-15 16:03:39.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.1845
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.1771

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000

Gjam

data<-sim_data$CompDenseSp10

#fit_gjam(data,5000,500,"./gjam_models/gjam10cmpd.rda",interact=comp_inter[[11]])

load_gjam(data,name="./gjam_models/gjam10cmpd.rda", interact= (-1)*comp_inter[[11]])
## 
## Sensitivity by predictor variables f:
##      Estimate    SE CI_025 CI_975
## env      88.6 22.30   52.1  138.0
## env2     40.9  9.14   24.7   60.8
## 
## Coefficient matrix B:
##             sp01   sp02   sp03   sp04   sp05   sp06   sp07   sp08  sp09
## intercept -0.123 -1.360  0.287  1.150  0.552  1.230 -0.933  0.459 0.629
## env       -2.520 -0.489 -1.150 -1.210 -0.255  0.376  0.868  1.790 2.790
## env2       0.141  0.308 -0.995 -0.557 -0.577 -0.833 -0.830 -0.579 0.276
## RMSPE      0.297  0.341  0.389  0.327  0.493  0.385  0.450  0.318 0.301
##             sp10
## intercept -0.519
## env        1.370
## env2       0.513
## RMSPE      0.359
## 
## Coefficient matrix B:
##                Estimate     SE  CI_025  CI_975 sig95
## sp01_intercept   -0.123 0.0954 -0.2940  0.0708      
## sp01_env         -2.520 0.1750 -2.8800 -2.2100     *
## sp01_env2         0.141 0.1240 -0.0623  0.3830      
## sp02_intercept   -1.360 0.0758 -1.5100 -1.2100     *
## sp02_env         -0.489 0.0590 -0.6020 -0.3750     *
## sp02_env2         0.308 0.0741  0.1760  0.4580     *
## sp03_intercept    0.287 0.0798  0.1340  0.4330     *
## sp03_env         -1.150 0.0816 -1.3100 -1.0000     *
## sp03_env2        -0.995 0.0941 -1.1600 -0.8190     *
## sp04_intercept    1.150 0.0925  0.9740  1.3300     *
## sp04_env         -1.210 0.1000 -1.4000 -1.0200     *
## sp04_env2        -0.557 0.0767 -0.7390 -0.4220     *
## sp05_intercept    0.552 0.0611  0.4330  0.6720     *
## sp05_env         -0.255 0.0535 -0.3610 -0.1530     *
## sp05_env2        -0.577 0.0637 -0.7010 -0.4560     *
## sp06_intercept    1.230 0.0730  1.0900  1.3800     *
## sp06_env          0.376 0.0650  0.2530  0.4990     *
## sp06_env2        -0.833 0.0821 -0.9990 -0.6810     *
## sp07_intercept   -0.933 0.0741 -1.0800 -0.7860     *
## sp07_env          0.868 0.0745  0.7100  1.0000     *
## sp07_env2        -0.830 0.0831 -1.0100 -0.6750     *
## sp08_intercept    0.459 0.0739  0.3310  0.6170     *
## sp08_env          1.790 0.1070  1.5900  2.0100     *
## sp08_env2        -0.579 0.0596 -0.6980 -0.4630     *
## sp09_intercept    0.629 0.0599  0.5130  0.7470     *
## sp09_env          2.790 0.1360  2.5400  3.0700     *
## sp09_env2         0.276 0.0597  0.1570  0.3880     *
## sp10_intercept   -0.519 0.0622 -0.6470 -0.4040     *
## sp10_env          1.370 0.1000  1.1800  1.5600     *
## sp10_env2         0.513 0.0589  0.3890  0.6240     *
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X and W:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
##  Design Table
##      env env2
## VIF    1    1
## env2   0   NA
## 
##  Sample contains n = 500 observations on S = 10 response variables, and 2 predictors.  Data types (typeNames) include PA. There are 0 missing values in X and 0 missing values in Y. The RMSPE is 0.371, and the DIC is 83249.  Computation involved 5000 Gibbs steps, with a burnin of 500.

#gjfd5<-load_object("./gjam_models/gjam10cmpd.rda")

#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))

HMSC

data<-sim_data$CompDenseSp10
hm_mod<-fit_hmsc(data,"Load",nsamples=5000, nchains=2,name="./HMmodels/hm10cmpd.rda" )
hm_conv(hm_mod)

hm_inter(hm_mod,nsamples=5000, nchains=2, interact =  (-1)*comp_inter[[11]])

Environmental filtering + Competition dense 20 species

## Summary for model '/tmp/RtmpKix4lZ/file40e96843124a'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 2061.424 minutes at time 2019-04-16 03:39:17.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.9065
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.6944

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000

Gjam

data<-sim_data$CompDenseSp20

#fit_gjam(data,10000,1000,"./gjam_models/gjam20cmpd.rda",interact= (-1)*comp_inter[[12]])

load_gjam(data,name="./gjam_models/gjam20cmpd.rda", interact= (-1)*comp_inter[[12]])
## 
## Sensitivity by predictor variables f:
##      Estimate   SE CI_025 CI_975
## env       318 54.6  218.0    432
## env2      116 20.8   78.2    160
## 
## Coefficient matrix B:
##              sp01   sp02   sp03   sp04   sp05   sp06    sp07   sp08   sp09
## intercept -0.0358 -1.130 -0.096  0.511 -0.184  0.893  0.0506  1.420  1.120
## env       -2.1300 -1.070 -1.910 -2.900 -1.300 -1.630 -0.5920 -0.832 -0.395
## env2       0.4330 -0.470 -0.696 -0.240 -0.649 -0.591 -0.3470 -0.637 -0.657
## RMSPE      0.3040  0.453  0.348  0.291  0.438  0.316  0.5230  0.343  0.394
##              sp10    sp11   sp12   sp13   sp14   sp15   sp16   sp17  sp18
## intercept  0.6690  0.0535  1.320  0.851  0.300  0.884 1.1300 -0.200 0.501
## env       -0.0441 -0.0424  0.831  0.931  0.686  1.690 2.5100  1.460 2.410
## env2      -0.4730 -0.2160 -0.519 -0.597 -0.967 -0.316 0.0994 -0.107 0.467
## RMSPE      0.4950  0.5700  0.355  0.388  0.445  0.331 0.3100  0.411 0.326
##              sp19   sp20
## intercept -0.4190 -2.040
## env        2.2000  0.776
## env2      -0.0511 -1.040
## RMSPE      0.3290  0.323
## 
## Coefficient matrix B:
##                Estimate     SE   CI_025  CI_975 sig95
## sp01_intercept  -0.0358 0.0820 -0.19400  0.1180      
## sp01_env        -2.1300 0.1020 -2.33000 -1.9400     *
## sp01_env2        0.4330 0.0873  0.26500  0.5990     *
## sp02_intercept  -1.1300 0.0847 -1.29000 -0.9620     *
## sp02_env        -1.0700 0.0834 -1.24000 -0.9120     *
## sp02_env2       -0.4700 0.1200 -0.69000 -0.2480     *
## sp03_intercept  -0.0960 0.0872 -0.26300  0.0815      
## sp03_env        -1.9100 0.1060 -2.13000 -1.7100     *
## sp03_env2       -0.6960 0.1170 -0.89800 -0.4690     *
## sp04_intercept   0.5110 0.0936  0.32900  0.6790     *
## sp04_env        -2.9000 0.1400 -3.18000 -2.6300     *
## sp04_env2       -0.2400 0.2230 -0.62100  0.0949      
## sp05_intercept  -0.1840 0.0748 -0.33600 -0.0353     *
## sp05_env        -1.3000 0.0866 -1.49000 -1.1400     *
## sp05_env2       -0.6490 0.1120 -0.89000 -0.4670     *
## sp06_intercept   0.8930 0.0924  0.70000  1.0600     *
## sp06_env        -1.6300 0.1040 -1.84000 -1.4400     *
## sp06_env2       -0.5910 0.0673 -0.73300 -0.4680     *
## sp07_intercept   0.0506 0.0585 -0.06590  0.1630      
## sp07_env        -0.5920 0.0617 -0.71400 -0.4700     *
## sp07_env2       -0.3470 0.0604 -0.46100 -0.2210     *
## sp08_intercept   1.4200 0.0957  1.23000  1.5900     *
## sp08_env        -0.8320 0.0847 -0.99900 -0.6660     *
## sp08_env2       -0.6370 0.0756 -0.78600 -0.4930     *
## sp09_intercept   1.1200 0.0746  0.97500  1.2700     *
## sp09_env        -0.3950 0.0564 -0.50700 -0.2870     *
## sp09_env2       -0.6570 0.0945 -0.84500 -0.4880     *
## sp10_intercept   0.6690 0.0746  0.51100  0.8000     *
## sp10_env        -0.0441 0.0569 -0.15200  0.0702      
## sp10_env2       -0.4730 0.0697 -0.61300 -0.3430     *
## sp11_intercept   0.0535 0.0655 -0.07700  0.1790      
## sp11_env        -0.0424 0.0578 -0.15800  0.0683      
## sp11_env2       -0.2160 0.0679 -0.34500 -0.0781     *
## sp12_intercept   1.3200 0.0891  1.14000  1.4800     *
## sp12_env         0.8310 0.0726  0.69500  0.9750     *
## sp12_env2       -0.5190 0.0947 -0.68700 -0.3440     *
## sp13_intercept   0.8510 0.0630  0.72600  0.9720     *
## sp13_env         0.9310 0.0732  0.77700  1.0700     *
## sp13_env2       -0.5970 0.0705 -0.73700 -0.4640     *
## sp14_intercept   0.3000 0.0604  0.19300  0.4350     *
## sp14_env         0.6860 0.0675  0.55500  0.8180     *
## sp14_env2       -0.9670 0.0751 -1.11000 -0.8230     *
## sp15_intercept   0.8840 0.0657  0.75500  1.0100     *
## sp15_env         1.6900 0.1250  1.48000  1.9600     *
## sp15_env2       -0.3160 0.0772 -0.46100 -0.1650     *
## sp16_intercept   1.1300 0.0883  0.97200  1.3200     *
## sp16_env         2.5100 0.1230  2.29000  2.7700     *
## sp16_env2        0.0994 0.0559 -0.00627  0.2120      
## sp17_intercept  -0.2000 0.0714 -0.33200 -0.0553     *
## sp17_env         1.4600 0.0885  1.29000  1.6300     *
## sp17_env2       -0.1070 0.0772 -0.26400  0.0315      
## sp18_intercept   0.5010 0.0802  0.35400  0.6640     *
## sp18_env         2.4100 0.1350  2.14000  2.6600     *
## sp18_env2        0.4670 0.0959  0.29100  0.6590     *
## sp19_intercept  -0.4190 0.0794 -0.57000 -0.2740     *
## sp19_env         2.2000 0.1060  2.00000  2.4100     *
## sp19_env2       -0.0511 0.0801 -0.20500  0.1060      
## sp20_intercept  -2.0400 0.1240 -2.28000 -1.8100     *
## sp20_env         0.7760 0.0763  0.63800  0.9350     *
## sp20_env2       -1.0400 0.0822 -1.20000 -0.8790     *
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X and W:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
##  Design Table
##      env env2
## VIF    1    1
## env2   0   NA
## 
##  Sample contains n = 500 observations on S = 20 response variables, and 2 predictors.  Data types (typeNames) include PA. There are 0 missing values in X and 0 missing values in Y. The RMSPE is 0.392, and the DIC is 431994.  Computation involved 10000 Gibbs steps, with a burnin of 1000.

#gjfd5<-load_object("./gjam_models/gjam20cmpd.rda")

#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))

HMSC

data<-sim_data$CompDenseSp20

hm_mod<-fit_hmsc(data,"Load",nsamples=5000, nchains=2,name="./HMmodels/hm20cmpd.rda" )
hm_conv(hm_mod)

hm_inter(hm_mod, nsamples=5000, nchains=2,interact =  (-1)*comp_inter[[12]])

Environmental filtering + Competition sparse 5 species

## Summary for model '/tmp/RtmpKix4lZ/file40e97e225117'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 29.251 minutes at time 2019-04-17 14:00:46.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.1876
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.2038

## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000

Gjam

data<-sim_data$CompSparseSp5


#fit_gjam(data,5000,2000,"./gjam_models/gjam5cmps.rda",interact= (-1)*comp_inter[[13]])

load_gjam(data,name="./gjam_models/gjam5cmps.rda", interact= (-1)*comp_inter[[13]])
## 
## Sensitivity by predictor variables f:
##      Estimate   SE CI_025 CI_975
## env      39.5 14.7   16.3   77.2
## env2     48.7 14.0   24.0   78.8
## 
## Coefficient matrix B:
##              sp01   sp02   sp03  sp04   sp05
## intercept -1.4600 -0.221  1.440 2.070 -0.115
## env       -0.9350 -1.160 -0.284 3.070  2.070
## env2       0.0755 -1.220 -0.798 0.521  0.528
## RMSPE      0.3580  0.423  0.357 0.336  0.329
## 
## Coefficient matrix B:
##                Estimate     SE  CI_025 CI_975 sig95
## sp01_intercept  -1.4600 0.0769 -1.6200 -1.320     *
## sp01_env        -0.9350 0.0824 -1.1100 -0.781     *
## sp01_env2        0.0755 0.0610 -0.0379  0.197      
## sp02_intercept  -0.2210 0.0475 -0.3150 -0.125     *
## sp02_env        -1.1600 0.0737 -1.3100 -1.020     *
## sp02_env2       -1.2200 0.0753 -1.3700 -1.080     *
## sp03_intercept   1.4400 0.0860  1.2700  1.610     *
## sp03_env        -0.2840 0.0584 -0.3930 -0.165     *
## sp03_env2       -0.7980 0.0661 -0.9300 -0.669     *
## sp04_intercept   2.0700 0.3490  1.5300  2.690     *
## sp04_env         3.0700 0.3970  2.4300  3.740     *
## sp04_env2        0.5210 0.1090  0.3170  0.736     *
## sp05_intercept  -0.1150 0.1890 -0.4600  0.291      
## sp05_env         2.0700 0.1450  1.8000  2.350     *
## sp05_env2        0.5280 0.1950  0.1750  0.946     *
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X and W:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
##  Design Table
##      env env2
## VIF    1    1
## env2   0   NA
## 
##  Sample contains n = 500 observations on S = 5 response variables, and 2 predictors.  Data types (typeNames) include PA. There are 0 missing values in X and 0 missing values in Y. The RMSPE is 0.362, and the DIC is 45260.  Computation involved 5000 Gibbs steps, with a burnin of 2000.

#gjfd5<-load_object("./gjam_models/gjam5cmps.rda")

#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam5cmps.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))

HMSC

data<-sim_data$CompSparseSp5

hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hm5cmps.rda" )
hm_conv(hm_mod)

hm_inter(hm_mod, nsamples=2000, nchains=2,interact =  (-1)*comp_inter[[13]])

Environmental filtering + Competition sparse 10 species

Gjam

data<-sim_data$CompSparseSp10
#fit_gjam(data,2000,500,"./gjam_models/gjam10cmps.rda",interact= (-1)*comp_inter[[14]])
load_gjam(data,name="./gjam_models/gjam10cmps.rda", interact= (-1)*comp_inter[[14]])
## 
## Sensitivity by predictor variables f:
##      Estimate    SE CI_025 CI_975
## env        82 20.10   44.5  124.0
## env2       45  9.55   31.6   67.9
## 
## Coefficient matrix B:
##             sp01    sp02   sp03   sp04   sp05   sp06    sp07    sp08  sp09
## intercept -1.220  0.0536  0.645  0.406  1.420  1.110  0.0665 -0.1640 0.887
## env       -1.670 -1.9900 -1.910 -0.238 -0.354  0.294  0.3120  1.0900 3.340
## env2       0.508  0.0778 -0.257 -0.998 -1.060 -0.581 -0.9340 -0.0165 0.105
## RMSPE      0.309  0.3470  0.341  0.436  0.354  0.393  0.4810  0.4600 0.277
##             sp10
## intercept -0.023
## env        2.390
## env2       0.548
## RMSPE      0.315
## 
## Coefficient matrix B:
##                Estimate     SE    CI_025  CI_975 sig95
## sp01_intercept  -1.2200 0.0857 -1.380000 -1.0400     *
## sp01_env        -1.6700 0.1090 -1.870000 -1.4500     *
## sp01_env2        0.5080 0.0571  0.395000  0.6200     *
## sp02_intercept   0.0536 0.0694 -0.067400  0.1920      
## sp02_env        -1.9900 0.1090 -2.200000 -1.7700     *
## sp02_env2        0.0778 0.0677 -0.053200  0.2100      
## sp03_intercept   0.6450 0.0544  0.538000  0.7490     *
## sp03_env        -1.9100 0.1290 -2.180000 -1.6800     *
## sp03_env2       -0.2570 0.0484 -0.350000 -0.1630     *
## sp04_intercept   0.4060 0.0555  0.295000  0.5140     *
## sp04_env        -0.2380 0.0728 -0.373000 -0.0994     *
## sp04_env2       -0.9980 0.0661 -1.130000 -0.8770     *
## sp05_intercept   1.4200 0.0719  1.280000  1.5600     *
## sp05_env        -0.3540 0.0592 -0.477000 -0.2440     *
## sp05_env2       -1.0600 0.0625 -1.180000 -0.9430     *
## sp06_intercept   1.1100 0.0611  0.981000  1.2200     *
## sp06_env         0.2940 0.0505  0.194000  0.3920     *
## sp06_env2       -0.5810 0.0517 -0.681000 -0.4820     *
## sp07_intercept   0.0665 0.0581 -0.049200  0.1770      
## sp07_env         0.3120 0.0594  0.193000  0.4190     *
## sp07_env2       -0.9340 0.0652 -1.060000 -0.8040     *
## sp08_intercept  -0.1640 0.0490 -0.258000 -0.0667     *
## sp08_env         1.0900 0.0732  0.944000  1.2200     *
## sp08_env2       -0.0165 0.0553 -0.129000  0.0802      
## sp09_intercept   0.8870 0.1040  0.705000  1.0800     *
## sp09_env         3.3400 0.1690  3.030000  3.6900     *
## sp09_env2        0.1050 0.0542  0.000574  0.2170     *
## sp10_intercept  -0.0230 0.0613 -0.134000  0.1000      
## sp10_env         2.3900 0.1260  2.160000  2.6500     *
## sp10_env2        0.5480 0.0675  0.405000  0.6680     *
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X and W:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
##  Design Table
##      env env2
## VIF    1    1
## env2   0   NA
## 
##  Sample contains n = 500 observations on S = 10 response variables, and 2 predictors.  Data types (typeNames) include PA. There are 0 missing values in X and 0 missing values in Y. The RMSPE is 0.377, and the DIC is 96021.  Computation involved 2000 Gibbs steps, with a burnin of 500.

#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam10cmps.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
data<-sim_data$CompSparseSp10
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hm10cmps.rda" )
hm_conv(hm_mod)

hm_inter(hm_mod, nsamples=2000, nchains=2,interact =  (-1)*comp_inter[[14]])

Environmental filtering + Competition sparse 20 species

Gjam

data<-sim_data$CompSparseSp20
#fit_gjam(data,2000,500,"./gjam_models/gjam20cmps.rda",interact= (-1)*comp_inter[[15]])
load_gjam(data,name="./gjam_models/gjam20cmps.rda", interact= (-1)*comp_inter[[15]])
## 
## Sensitivity by predictor variables f:
##      Estimate   SE CI_025 CI_975
## env       288 84.9    156    457
## env2      128 21.3     94    175
## 
## Coefficient matrix B:
##             sp01   sp02   sp03   sp04    sp05   sp06   sp07   sp08   sp09
## intercept -0.366 -0.677 -0.599  0.151  0.2260  0.817  0.510  0.979  1.070
## env       -1.990 -1.940 -2.400 -1.920 -1.7600 -1.760 -0.949 -1.030 -0.562
## env2       0.777 -0.456 -0.564 -0.531  0.0776 -0.247 -0.489 -0.535 -0.856
## RMSPE      0.316  0.361  0.321  0.329  0.3680  0.346  0.423  0.364  0.386
##             sp10    sp11   sp12   sp13   sp14   sp15    sp16   sp17
## intercept  1.020  0.8200  1.180  1.050  1.090  0.757  0.7880  0.479
## env       -0.291  0.0736  0.465  0.842  1.330  1.520  1.8300  2.380
## env2      -0.928 -0.6110 -0.634 -0.547 -0.554 -0.353 -0.0702 -0.104
## RMSPE      0.421  0.4590  0.383  0.383  0.330  0.356  0.3450  0.317
##              sp18    sp19   sp20
## intercept  0.0281 -0.2080 -0.808
## env        1.9000  2.0400  2.280
## env2      -0.0839  0.0373 -0.101
## RMSPE      0.3490  0.3430  0.328
## 
## Coefficient matrix B:
##                Estimate     SE  CI_025   CI_975 sig95
## sp01_intercept  -0.3660 0.0984 -0.5770 -0.20900     *
## sp01_env        -1.9900 0.1510 -2.2600 -1.71000     *
## sp01_env2        0.7770 0.0681  0.6550  0.91500     *
## sp02_intercept  -0.6770 0.0865 -0.8320 -0.49100     *
## sp02_env        -1.9400 0.0850 -2.1100 -1.78000     *
## sp02_env2       -0.4560 0.0532 -0.5570 -0.35000     *
## sp03_intercept  -0.5990 0.0619 -0.7310 -0.48600     *
## sp03_env        -2.4000 0.1150 -2.6200 -2.18000     *
## sp03_env2       -0.5640 0.0586 -0.6800 -0.45600     *
## sp04_intercept   0.1510 0.0661  0.0233  0.27700     *
## sp04_env        -1.9200 0.1190 -2.1200 -1.67000     *
## sp04_env2       -0.5310 0.0597 -0.6460 -0.41500     *
## sp05_intercept   0.2260 0.0559  0.1200  0.33300     *
## sp05_env        -1.7600 0.0851 -1.9400 -1.60000     *
## sp05_env2        0.0776 0.0593 -0.0411  0.19100      
## sp06_intercept   0.8170 0.0581  0.6950  0.92700     *
## sp06_env        -1.7600 0.0971 -1.9500 -1.58000     *
## sp06_env2       -0.2470 0.0566 -0.3600 -0.14000     *
## sp07_intercept   0.5100 0.0540  0.4030  0.61700     *
## sp07_env        -0.9490 0.0571 -1.0600 -0.83000     *
## sp07_env2       -0.4890 0.0590 -0.6070 -0.37800     *
## sp08_intercept   0.9790 0.0615  0.8600  1.10000     *
## sp08_env        -1.0300 0.0912 -1.1900 -0.85900     *
## sp08_env2       -0.5350 0.0543 -0.6450 -0.42800     *
## sp09_intercept   1.0700 0.0667  0.9430  1.21000     *
## sp09_env        -0.5620 0.0685 -0.6940 -0.43000     *
## sp09_env2       -0.8560 0.0573 -0.9660 -0.74400     *
## sp10_intercept   1.0200 0.0626  0.9000  1.15000     *
## sp10_env        -0.2910 0.0502 -0.3850 -0.19200     *
## sp10_env2       -0.9280 0.0578 -1.0400 -0.81900     *
## sp11_intercept   0.8200 0.0553  0.7100  0.93000     *
## sp11_env         0.0736 0.0500 -0.0251  0.16800      
## sp11_env2       -0.6110 0.0570 -0.7290 -0.50000     *
## sp12_intercept   1.1800 0.0596  1.0600  1.30000     *
## sp12_env         0.4650 0.0592  0.3480  0.58300     *
## sp12_env2       -0.6340 0.0529 -0.7400 -0.53500     *
## sp13_intercept   1.0500 0.0582  0.9290  1.16000     *
## sp13_env         0.8420 0.0529  0.7340  0.93700     *
## sp13_env2       -0.5470 0.0493 -0.6370 -0.44800     *
## sp14_intercept   1.0900 0.0690  0.9610  1.22000     *
## sp14_env         1.3300 0.0651  1.2000  1.46000     *
## sp14_env2       -0.5540 0.0574 -0.6730 -0.44200     *
## sp15_intercept   0.7570 0.0580  0.6440  0.86900     *
## sp15_env         1.5200 0.0708  1.3800  1.66000     *
## sp15_env2       -0.3530 0.0507 -0.4540 -0.25600     *
## sp16_intercept   0.7880 0.0591  0.6680  0.89500     *
## sp16_env         1.8300 0.1020  1.6300  2.02000     *
## sp16_env2       -0.0702 0.0457 -0.1550  0.02360      
## sp17_intercept   0.4790 0.0492  0.3830  0.56900     *
## sp17_env         2.3800 0.0979  2.1900  2.58000     *
## sp17_env2       -0.1040 0.0494 -0.2020 -0.00328     *
## sp18_intercept   0.0281 0.0502 -0.0737  0.12500      
## sp18_env         1.9000 0.1020  1.7100  2.10000     *
## sp18_env2       -0.0839 0.0534 -0.1840  0.02740      
## sp19_intercept  -0.2080 0.0626 -0.3440 -0.09310     *
## sp19_env         2.0400 0.1380  1.7600  2.29000     *
## sp19_env2        0.0373 0.0524 -0.0618  0.13600      
## sp20_intercept  -0.8080 0.0647 -0.9280 -0.67900     *
## sp20_env         2.2800 0.1060  2.0700  2.48000     *
## sp20_env2       -0.1010 0.0547 -0.2090  0.00857      
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X and W:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
##  Design Table
##      env env2
## VIF    1    1
## env2   0   NA
## 
##  Sample contains n = 500 observations on S = 20 response variables, and 2 predictors.  Data types (typeNames) include PA. There are 0 missing values in X and 0 missing values in Y. The RMSPE is 0.363, and the DIC is 340096.  Computation involved 2000 Gibbs steps, with a burnin of 500.

#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam20cmps.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))

HMSC

data<-sim_data$CompSparseSp20

hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hm20cmps.rda" )
hm_conv(hm_mod)

hm_inter(hm_mod, nsamples=2000, nchains=2,interact =  (-1)*comp_inter[[15]])

Environmental filtering + Competition +Facilitation 5 dense species

Environmental filtering + Competition +Facilitation 10 dense species

Environmental filtering + Competition +Facilitation 20 dense species

Graph from paper

Visualization for true interactions